function [mask, weight] = mesh2mask(node, face, xi, yi, hf, varargin)
%
% [mask weight]=mesh2mask(node,face,Nxy)
%   or
% [mask weight]=mesh2mask(node,face,[Nx,Ny])
%   or
% [mask weight]=mesh2mask(node,face,xi,yi,hf)
%
% fast rasterization of a 2D mesh to an image with triangle index labels
%
% author: Qianqian Fang <q.fang at neu.edu>
% date for initial version: July 18,2013
%
% input:
%      node: node coordinates, dimension N by 2 or N by 3 array
%      face: a triangle surface, N by 3 or N by 4 array
%      Nx,Ny,Nxy: output image in x/y dimensions, or both
%      xi,yi: linear vectors for the output pixel center positions in x/y
%      hf: (optional) the handle of a pre-created figure window, for faster
%          rendering
%
% output:
%      mask: a 2D image, the value of each pixel is the index of the
%            enclosing triangle, if the pixel is outside of the mesh, NaN
%      weight: (optional) a 3 by Nx by Ny array, where Nx/Ny are the dimensions for
%            the mask
%
% note: This function only works in MATLAB when the DISPLAY is not
%       disabled. The maximum size of the mask output is limited by the
%       screen size.
%
% example:
%
%   [no,fc]=meshgrid6(0:5,0:5);
%   [mask weight]=mesh2mask(no,fc,-1:0.1:5,0:0.1:5);
%   imagesc(mask);
%
% -- this function is part of iso2mesh toolbox (http://iso2mesh.sf.net)
%

opt = varargin2struct(varargin{:});

if (nargin == 3 && length(xi) == 1 && xi > 0)
    mn = min(node);
    mx = max(node);
    df = (mx(1:2) - mn(1:2)) / xi;
elseif (nargin == 3 && length(xi) == 2 && all(xi > 0))
    mn = min(node);
    mx = max(node);
    df = (mx(1:2) - mn(1:2)) ./ xi;
elseif (nargin == 4 || nargin == 5)
    mx = [max(xi) max(yi)];
    mn = [min(xi) min(yi)];
    df = [min(diff(xi(:))) min(diff(yi(:)))];
else
    error('you must give at least xi input');
end
if (size(node, 2) <= 1 || size(face, 2) <= 2)
    error('node must have 2 or 3 columns; face can not have less than 2 columns');
end

if (nargin < 5)
    hf = figure('visible', 'on');
else
    clf(hf);
end

if (jsonopt('edge', 0, opt))
    try
        set(gcf, 'GraphicsSmoothing', 'off');
    catch
    end
end

patch('Vertices', node, 'Faces', face, 'linestyle', 'none', 'FaceColor', 'flat', ...
      'FaceVertexCData', (1:size(face, 1))', 'CDataMapping', 'scaled');
set(gca, 'Position', [0 0 1 1]);
cm = jet(size(face, 1));
colormap(cm);
axis off;
set(gca, 'xlim', [mn(1) mx(1)]);
set(gca, 'ylim', [mn(2) mx(2)]);
set(gca, 'clim', [1 size(face, 1)]);

output_size = round((mx(1:2) - mn(1:2)) ./ df); % Size in pixels

if (isoctavemesh || isempty(getenv('DISPLAY')))
    resolution = 300; % Resolution in DPI
    set(gcf, 'PaperPositionMode', 'manual');
    set(gcf, 'paperunits', 'inches', 'paperposition', [0 0 output_size / resolution]);
    deletemeshfile(mwpath('post_mesh2mask.png'));
    print(mwpath('post_mesh2mask.png'), '-dpng', ['-r' num2str(resolution)]);
    mask = imread(mwpath('post_mesh2mask.png'));
else
    pos = get(hf, 'position');
    pos(3:4) = max(pos(3:4), output_size + 20);
    set(hf, 'position', pos);
    set(gca, 'Units', 'pixels', 'position', [1, 1, output_size(1), output_size(2)]);
    mask = getframe(gca);
    if (any(size(mask.cdata) < [output_size([2 1]) 3]))
        error('the requested rasterization grid is larger than the screen resolution');
    end
    mask = mask.cdata((end - output_size(2) + 1):end, 1:output_size(1), :);
end
if (nargin < 5)
    close(hf);
end
mask = int32(reshape(mask, [size(mask, 1) * size(mask, 2) size(mask, 3)]));
[isfound, locb] = ismember(mask, floor(cm * 255), 'rows');
locb(isfound == 0) = nan;

mask = rot90(reshape(locb, output_size([2 1]))');

if (nargout >= 2)
    xi = mn(1) + df(1) / 2:df(1):mx(1);
    yi = mn(2) + df(2) / 2:df(2):mx(2);
    weight = barycentricgrid(node, face, xi, yi, mask);
    if (size(face, 2) >= 4)
        badidx = find(weight(1, :, :) < 0 | weight(2, :, :) < 0 | weight(3, :, :) < 0);
        badidx = badidx(face(mask(badidx), 3) ~= face(mask(badidx), 4));
        weight2 = barycentricgrid(node, face(:, [1 3 4]), xi, yi, mask);
        weight(:, badidx) = 0;
        weight([1 3 4], badidx) = weight2(:, badidx);
    end
end

function weight = barycentricgrid(node, face, xi, yi, mask)
[xx, yy] = meshgrid(xi, yi);
idx = find(~isnan(mask));
eid = mask(idx);
t1 = node(face(:, 1), :);
t2 = node(face(:, 2), :);
t3 = node(face(:, 3), :);
tt = (t2(:, 2) - t3(:, 2)) .* (t1(:, 1) - t3(:, 1)) + (t3(:, 1) - t2(:, 1)) .* (t1(:, 2) - t3(:, 2));
w(:, 1) = (t2(eid, 2) - t3(eid, 2)) .* (xx(idx) - t3(eid, 1)) + (t3(eid, 1) - t2(eid, 1)) .* (yy(idx) - t3(eid, 2));
w(:, 2) = (t3(eid, 2) - t1(eid, 2)) .* (xx(idx) - t3(eid, 1)) + (t1(eid, 1) - t3(eid, 1)) .* (yy(idx) - t3(eid, 2));
w(:, 1) = w(:, 1) ./ tt(eid);
w(:, 2) = w(:, 2) ./ tt(eid);
w(:, 3) = 1 - w(:, 1) - w(:, 2);
weight = zeros(3, size(mask, 1), size(mask, 2));
ww = zeros(size(mask));
ww(idx) = w(:, 1);
weight(1, :, :) = ww;
ww(idx) = w(:, 2);
weight(2, :, :) = ww;
ww(idx) = w(:, 3);
weight(3, :, :) = ww;
